function [gamma,gamma_x,gamma_y,gamma_z] = Threaded3DGamma(RefDose,EvalDose,VoxelLength,DoseDiff,DTA,nThreads,interp)

% Threaded3DGamma
% Based on BrutalGamma2d by Michael Jensen.
% Splits the task of Gamma calculation up to multiple threads.
% Computes Gamma for 3D volumes.

% Inputs:
%   RefDose - Reference dose matrix (the points for which you get a gamma
%             back)
%   EvalDose - Evaluation dose matrix (the one that will be interpolated)
%   VoxelLength - the length of the cubic voxel in mm
%   DoseDiff - the dose difference criterion in % of maximum dose 
%              eg 3 for 3%
%   DTA - the distance to agreement criterion in mm
%         eg 3  for 3mm
%   nThreads - Number of threads to spawn.
%   interp - % of DTA to interpolate evaluation dose matrix down to. Optional
%           eg 0.10 for 1/10th the DTA criterion.

% Outputs:
%   gamma - 3D matrix with gamma values
%   gamma_x - 3D matrix with displacement values in X axis
%   gamma_y - 3D matrix with displacement values in y axis
%   gamma_z - 3D matrix with displacement values in z axis
%		Displacement values are only for gamma values over 1, where there is a displacement.

% Michael D. Jensen
% The University of Western Ontario
% December 2010
%Modified by:
% Jeff Craig
% July 2011


% Copyright 2011, Jeff Craig and Michael Jensen.
%  
% Threaded 3D Gamma is distributed under the terms of the Lesser GNU Public License. 
% 
%     This version of Threaded 3D Gamma is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
% Threaded 3D Gamma is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with Threaded 3D Gamma.  If not, see <http://www.gnu.org/licenses/>.


if nargin < 7
	interp = 1;
end


% Method used for 2D interpolation
interpMethod = 'Spline';

% Compute the dose difference criterion as percentage of maximum dose
dDose = DoseDiff*0.01*max(RefDose(:)); 

% Reference Dose and Evaluation Dose must be the same size
if size(RefDose)~=size(EvalDose)
    error('Dose Matrices are not the same size.')
end

% Reference Dose and Evaluation Dose must be 3D
% if ndims(RefDose)~=3
%     error('Dose Matrices are not 3D.')
% end

% Cannot spawn more threads then indices.
if length(RefDose(:))<nThreads
    error('To many threads.')
end

% Dose Matrix Dimensions
[nx ny nz] = size(EvalDose);

%make meshgrid for original resolution dose matrix
xruler = (1:nx)*VoxelLength;
yruler = (1:ny)*VoxelLength;
zruler = (1:nz)*VoxelLength;
[x y z] = ndgrid(xruler,yruler,zruler);

%determine if it is necessary to interpolate Eval dose
if interp == 1
    %Compute Gamma using mex function, for speed
    tic
    [gamma,gamma_x,gamma_y,gamma_z] = Threaded3DGammamex(RefDose,EvalDose,x,y,z,x,y,z,dDose,DTA,nThreads);
    toc
    
    figure();quiver3(x,y,z,gamma_x,gamma_y,gamma_z,0);
else
    %interpolate voxel size by fraction of the DTA criterion
    vx_interp = DTA * interp;
    
    %%% Interpolate Evaluation Dose to finer spacing
    %make interpolated mesh grid
    xxruler = vx_interp : vx_interp : xruler(end);
    yyruler = vx_interp : vx_interp : yruler(end);
    zzruler = vx_interp : vx_interp : zruler(end);
    [xx yy zz] = ndgrid(xxruler,yyruler,zzruler);
    
    tic
    %interpolate evaluation dose
    EvalDoseInterp = interpn(x,y,z,EvalDose,xx,yy,zz,interpMethod);
    toc
    
    disp('Finished Interpolation')
    
    %Compute Gamma using mex function, for speed
    tic
    [gamma,gamma_x,gamma_y,gamma_z] = Threaded3DGammamex(RefDose,EvalDoseInterp,x,y,z,xx,yy,zz,dDose,DTA,nThreads);
    toc
    
    figure();quiver3(x,y,z,gamma_x,gamma_y,gamma_z,0);
end

end
